$ontext
This code is for the following paper:
    Cai, Yongyang, William Brock, and Anastasios Xepapadeas (2022).
    Climate Change Impact on Economic Growth: Regional Climate Policy under Cooperation and Noncooperation.
    JAERE, forthcoming

Author: Yongyang Cai, The Ohio State University
Date: August 11, 2022
$offtext

parameter elasmu elasticity of marginal utility / 1.45 /
    nlag    / 10 /;

set t  Time periods (1 year per period)                   	/2015*2515/ ;
set tfirst(t);
tfirst(t)$(ord(t)=1) = YES;

alias(t,s);

set r /North, Tropic, South/;
set ii(r) /North, Tropic/;

parameters
** Preferences
    beta    discount factor per year   / .985  /

** Population and technology
    gama     Capital elasticity in production function        /.300    /
    dk       Depreciation rate on capital (per year)          /.100    /

** Emissions parameters
    eland0   Carbon emissions from land 2015 (GtC per year)     / 0.95    /
    deland   Decline rate of land emissions (per year)          / .115     /

** Carbon cycle
* Initial Conditions
    mat0   Initial Concentration in atmosphere 2015 (GtC)        /851   /
    mu0    Initial Concentration in upper strata 2015 (GtC)      /460   /
	ml0    Initial Concentration in lower strata 2015 (GtC)      /1740   /
    mateq  Equilibrium concentration atmosphere  (GtC)           /588     /
    mueq   Equilibrium concentration in upper strata (GtC)       /360    /
    mleq   Equilibrium concentration in lower strata (GtC)       /1720   /

** Climate model parameters
    fex0     2015 forcings of non-CO2 GHG (Wm-2)              / 0.5  /
    fex1     2100 forcings of non-CO2 GHG (Wm-2)              / 1  /

    fco22x   Forcings of equilibrium CO2 doubling (Wm-2)      /3.68    /

** Abatement cost
    expcost2  Exponent of control cost function               / 2.6  /
    gback     Initial cost decline backstop cost per year    / .005 /
	LIMMIU     Upper limit on control rate                     / 1.2    /
;


PARAMETERS
    sig0(ii)
    gsigma1(ii)  Initial growth of sigma (per year)                
    dsig(ii)     Decline rate of decarbonization (per year)        
    sigma(t,ii)      Carbon-equivalent-emissions output ratio
    cost1(t,ii)      Adjusted cost for backstop
    etree(t)      Emissions from deforestation
    forcoth(t)    Exogenous forcing for other greenhouse gases
    elasmu1
    pback(ii)     Cost of backstop 2005$ per tC 2010           
;

$gdxin carbonIntensity_param
$load sig0 gsigma1 dsig

sigma(t,ii) = sig0(ii)*exp(-gsigma1(ii)*(1-exp(-dsig(ii)*(ord(t)-1)))/dsig(ii));

* RICE data (in 2005 $), converted to 2010$ (2.23% per year in 2005-2010 according to the Bureau of Labor Statistics consumer price index)
pback('North') = 1.13 * 1.165;
pback('Tropic') = 1.44 * 1.165;

cost1(t,ii) = pback(ii)*exp(-gback*(ord(t)-1))*sigma(t,ii)/expcost2;
etree(t) = eland0*exp(-deland*(ord(t)-1));
forcoth(t) = fex0+ (1/85)*(fex1-fex0)*(ord(t)-1)$(ord(t) lt 86)+ (fex1-fex0)$(ord(t) ge 86);

elasmu1 = 1-elasmu;

PARAMETERS
* Flow paramaters
    b12      Carbon cycle transition matrix                      
	b11          Carbon cycle transition matrix
 	b21          Carbon cycle transition matrix
	b22          Carbon cycle transition matrix
	b23
	b32 
	b33

	k0(ii)       Initial capital value (trill 2010 USD) 
	q0(ii)		initial gross regional output (trill 2010 USD)   
	tocean0     Initial lower stratum temp change above the preindustrial level 
	tatm0(r)    Initial atmospheric temp change above the preindustrial level 
	a0(ii)       Initial level of total factor productivity       
	ga0(ii)      Initial growth rate for TFP per year             
	dela(ii)     Decline rate of TFP per year in the absence of climate damage     
	ga1(ii)	     temperature impact on TFP   
    ga2(ii)      temperature impact on TFP
    ga(t,ii)
    delaT(ii)     Decline rate of climate impact per year
    AL_noDam(t,ii)
    Pop(t,ii)

	xi1       Climate equation coefficient for upper level     
	xi2       Transfer coefficient upper to lower stratum     
	xi3       Transfer coefficient for lower level             
	B
	xi4
    xi5
    xi7
    xi8
;

$gdxin climate_anomaly_param
$load b12 b23 xi1 xi2 xi3 B xi4 xi5 xi7 xi8


set gd /ga0, dela/;
set gCC /ga1, ga2, delaT/;

$offlisting

Table Pop(t,ii) 
$ondelim
$include RegionalPop_SSP1_Annual.csv
$offdelim
;

Table paraTFPnoCC(gd,ii) 
$ondelim
$include paramTFP_noDam.csv
$offdelim
;

Table paraTFP_CC(gCC,ii) 
$ondelim
$include paramTFP_Dam_lag10.csv
$offdelim
;

$onlisting

ga0(ii) = paraTFPnoCC('ga0',ii);
dela(ii) = paraTFPnoCC('dela',ii);
ga1(ii) = paraTFP_CC('ga1',ii);
ga2(ii) = paraTFP_CC('ga2',ii);
delaT(ii) = paraTFP_CC('delaT',ii);

display dela, ga0, ga1, ga2, delaT;

b11 = 1 - b12;
b21 = b12*mateq/mueq;
b22 = 1 - b21 - b23;
b32 = b23*mueq/mleq;
b33 = 1 - b32;

* RICE data 
k0('North') = 100;
k0('Tropic') = 53;

* Use world bank data (2015 GDP) 
q0('North') = 54;
q0('Tropic') = 19;

* time periods in CMIP5 
set tt / 2006*2100 /;
set n sets of scenarios / RCP26, RCP45, RCP6, RCP85 /;

$offlisting

Table TA_RCP85(tt,r) 
$ondelim
$include RegionalSurfaceTempAnomaly_RCP85.csv
$offdelim
;

Table TLs(tt,n) 
$ondelim
$include OceanTempAnomaly.csv
$offdelim
;

$onlisting

* 2015 temperature anomaly (from OceanTempAnomaly.csv and RegionalSurfaceTempAnomaly_RCP85.csv)
tocean0 = TLs('2015','RCP85');
tatm0(r) = TA_RCP85('2015',r);

ga(t,ii)=ga0(ii)*exp(-dela(ii)*(ord(t)-1));

A0(ii) = sum(tfirst, q0(ii) / (k0(ii)**GAMA * Pop(tfirst,ii)**(1-GAMA)));
AL_noDam(tfirst,ii) = A0(ii);
loop((t,ii)$(ord(t)<card(t)),
  AL_noDam(t+1,ii) = AL_noDam(t,ii)/(1-ga(t,ii));
);

display AL_noDam;

*****************************************************************

positive VARIABLES
 K(t,ii)            Capital stock trillions US dollars
 MAT(t)          Carbon concentration in atmosphere GtC
 MU(t)           Carbon concentration in shallow oceans Gtc
 ML(t)
 Cratio(t,ii)            Consumption trillions US dollars
 miu(t,ii)          Emission control rate GHGs
 Y(t,ii)
 AL(t,ii)
 netY(t,ii)
 mitiCost(t,ii)
;

variables
 F(t)
 E(t,ii)            CO2-equivalent emissions GtC
 TATM(t,r)         Temperature of atmosphere in degrees C
 TOCEAN(t)       Temperature of lower oceans degrees C
 UTILITY
;

*****************************************************************

EQUATIONS

 KK(t,ii)            Capital balance equation
 MMAT(t)          Atmospheric concentration equation
 MMU(t)           Shallow ocean concentration
 MML(t)
 TATMEQ_Tropic(t)        Temperature-climate equation for atmosphere
 TATMEQ_North(t)        Temperature-climate equation for atmosphere
 TATMEQ_South(t)        Temperature-climate equation for atmosphere
 TOCEANEQ(t)      Temperature-climate equation for lower oceans

 TFPeq(t,ii)

 ForcingEq(t)
 EE(t,ii)            Emissions equation
 MitigationCostFun(t,ii)
 Output(t,ii)
 NetOutput(t,ii)

 UTIL             Objective function
;

** Equations of the model

KK(t,ii)$(ord(t)<card(t))..       K(t+1,ii) =e= (1-DK)*K(t,ii) + (1-Cratio(t,ii))*netY(t,ii);
MMAT(t)$(ord(t)<card(t))..        MAT(t+1)  =e= MAT(t)*b11+MU(t)*b21 + sum(ii, E(t,ii)) + ETREE(t);
MMU(t)$(ord(t)<card(t))..         MU(t+1)   =e= MAT(t)*b12+MU(t)*b22+ML(t)*b32;
MML(t)$(ord(t)<card(t))..         ML(t+1)   =e= MU(t)*b23+ML(t)*b33;

TATMEQ_North(t)$(ord(t)<card(t))..     	TATM(t+1,'North') =e= (1-B)*TATM(t,'North') - 
						xi2*(TATM(t,'North')-TOCEAN(t)) -
						xi4*(TATM(t,'North')-TATM(t,'Tropic')) + 
						(xi1+xi7)*F(t);
TATMEQ_Tropic(t)$(ord(t)<card(t))..     	TATM(t+1,'Tropic') =e= (1-B)*TATM(t,'Tropic') - 
						xi2*(TATM(t,'Tropic')-TOCEAN(t)) - 
						xi4/2*(TATM(t,'Tropic')-TATM(t,'North')) - 
                        xi5/2*(TATM(t,'Tropic')-TATM(t,'South')) +
						(xi1+xi8)*F(t);
TATMEQ_South(t)$(ord(t)<card(t))..     	TATM(t+1,'South') =e= (1-B)*TATM(t,'South') - 
						xi2*(TATM(t,'South')-TOCEAN(t)) -
						xi5*(TATM(t,'South')-TATM(t,'Tropic')) + 
						xi1*F(t);

TOCEANEQ(t)$(ord(t)<card(t))..   	TOCEAN(t+1) =e= TOCEAN(t) + sum(r, xi3*(TATM(t,r)-TOCEAN(t))) + xi3*(TATM(t,'Tropic')-TOCEAN(t));


ForcingEq(t)$(ord(t)<card(t))..   	F(t) =e= FCO22X*(log(MAT(t)/mateq)/log(2))+FORCOTH(t);

EE(t,ii)$(ord(t)<card(t))..          	E(t,ii) =e= sigma(t,ii) * (1-miu(t,ii))*Y(t,ii);

TFPeq(t,ii)$(ord(t)<card(t))..		AL(t,ii) =e= AL_noDam(t,ii) / ( 1+sum(s$(ord(s)<=ord(t) and ord(s)>=ord(t)-nlag),
                                        (ga1(ii)*(TATM(s,ii)-TATM0(ii))+ga2(ii)*sqr(TATM(s,ii)-TATM0(ii)))*delaT(ii)**(ord(t)-ord(s))) );

Output(t,ii)$(ord(t)<card(t))..		Y(t,ii) =e= (AL(t,ii)*Pop(t,ii)**(1-GAMA)*K(t,ii)**GAMA);
MitigationCostFun(t,ii)$(ord(t)<card(t)).. mitiCost(t,ii) =e= cost1(t,ii)*(miu(t,ii)**EXPcost2);
NetOutput(t,ii)$(ord(t)<card(t))..	netY(t,ii) =e= (1-mitiCost(t,ii))*Y(t,ii);

UTIL..             UTILITY =e= sum( (t,ii)$(ord(t)<card(t)), (Cratio(t,ii)*netY(t,ii)/Pop(t,ii))**elasmu1/elasmu1 * Pop(t,ii) * beta**(ord(t)-1) );

*******************************
**  Upper and Lower Bounds

K.lo(t,ii) = 0.01;
MAT.lo(t) = 10;
MU.lo(t) = 1;
Y.lo(t,ii) = 0.01;
netY.lo(t,ii) = 0.01;

Cratio.lo(t,ii) = 0.01;
Cratio.up(t,ii) = 1;
miu.lo(t,ii) = 0.001;
miu.up(t,ii) = LIMMIU;
miu.up(t,ii)$(ord(t)<150) = 1;

* Initial conditions
K.FX(tfirst,ii) = k0(ii);
MAT.FX(tfirst) = mat0;
MU.FX(tfirst) = mu0;
ML.fx(tfirst) = ml0;
TATM.FX(tfirst,r) = tatm0(r);
TOCEAN.FX(tfirst) = tocean0;
AL.fx(tfirst,ii) = A0(ii);

**********************************

** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = off;
option limrow = 0;
option limcol = 0;
option nlp = conopt;


model SocialProb / KK, MMAT, MMU, MML, TATMEQ_South, TATMEQ_Tropic, TATMEQ_North, TOCEANEQ,
  TFPeq, ForcingEq, MitigationCostFun, EE, Output, NetOutput, UTIL /;

parameters 
	nnn
	nsucc
	nnnmax / 20 /
	nsuccmax / 1 /
	opt_SCC(t,ii)
	opt_carbonTax(t,ii)
;

********************************************************************
* solve social planner's problem 
********************************************************************

miu.l(t,ii) = 0.99;
miu.l(t,'North')$(ord(t)<30) = 0.6 + (0.99-0.6)*(ord(t)-1)/29;
miu.l(t,'Tropic')$(ord(t)<70) = 0.16 + (0.99-0.16)*(ord(t)-1)/69;

mitiCost.l(t,ii) = cost1(t,ii)*(miu.l(t,ii)**EXPcost2);

Cratio.l(t,ii) = 0.75;

K.l(tfirst,ii) = k0(ii);
MAT.l(tfirst) = mat0;
MU.l(tfirst) = mu0;
ML.l(tfirst) = ml0;
TATM.l(tfirst,r) = tatm0(r);
TOCEAN.l(tfirst) = tocean0;
AL.l(tfirst,ii) = A0(ii);

loop(t$(ord(t)<card(t)),
    AL.l(t,ii) = AL_noDam(t,ii) * exp( -sum(s$(ord(s)<=ord(t) and ord(s)>=ord(t)-nlag),
        (ga1(ii)*(TATM.l(s,ii)-TATM0(ii))+ga2(ii)*sqr(TATM.l(s,ii)-TATM0(ii)))*delaT(ii)**(ord(t)-ord(s))) );

    Y.l(t,ii) = (AL.l(t,ii)*Pop(t,ii)**(1-GAMA)*K.l(t,ii)**GAMA);
    netY.l(t,ii) = (1-mitiCost.l(t,ii))*Y.l(t,ii);
    K.l(t+1,ii) = (1-DK)*K.l(t,ii) + Cratio.l(t,ii)*netY.l(t,ii);

    E.l(t,ii) = sigma(t,ii) * (1-miu.l(t,ii))*Y.l(t,ii);

    MAT.l(t+1)  = MAT.l(t)*b11+MU.l(t)*b21 + sum(ii, E.l(t,ii)) + ETREE(t);
    MU.l(t+1)   = MAT.l(t)*b12+MU.l(t)*b22+ML.l(t)*b32;
    ML.l(t+1)   = MU.l(t)*b23+ML.l(t)*b33;

    F.l(t) = FCO22X*(log(MAT.l(t)/mateq)/log(2))+FORCOTH(t);

    TATM.l(t+1,'North') = (1-B)*TATM.l(t,'North') - 
                            xi2*(TATM.l(t,'North')-TOCEAN.l(t)) -
                            xi4*(TATM.l(t,'North')-TATM.l(t,'Tropic')) + 
                            (xi1+xi7)*F.l(t);
    TATM.l(t+1,'Tropic') = (1-B)*TATM.l(t,'Tropic') - 
						xi2*(TATM.l(t,'Tropic')-TOCEAN.l(t)) - 
						xi4/2*(TATM.l(t,'Tropic')-TATM.l(t,'North')) - 
                        xi5/2*(TATM.l(t,'Tropic')-TATM.l(t,'South')) +
						(xi1+xi8)*F.l(t);
    TATM.l(t+1,'South') = (1-B)*TATM.l(t,'South') - 
						xi2*(TATM.l(t,'South')-TOCEAN.l(t)) -
						xi5*(TATM.l(t,'South')-TATM.l(t,'Tropic')) + 
						xi1*F.l(t);
						
    TOCEAN.l(t+1) = TOCEAN.l(t) + sum(r, xi3*(TATM.l(t,r)-TOCEAN.l(t))) + xi3*(TATM.l(t,'Tropic')-TOCEAN.l(t));
);
 
nnn = 0;
nsucc = 0;
option solprint = off;
Repeat(
    nnn = nnn+1;
    solve SocialProb maximizing UTILITY using nlp ;
    if((SocialProb.MODELSTAT eq 2 and SocialProb.SOLVESTAT eq 1),
      nsucc = nsucc+1;
    else
      nsucc = 0;
    );
until (nnn>=nnnmax or nsucc>=nsuccmax) );

ABORT$(nsucc < nsuccmax) "Model not normally completed";


*****************************

opt_SCC(t,ii)$(ord(t)<card(t))=-1*MMAT.m(t)*1000/(KK.m(t,ii)+0.00000000001)*12/44;
opt_carbonTax(t,ii)$(ord(t)<card(t))=(cost1(t,ii)*(expcost2*miu.l(t,ii)**(expcost2-1)/sigma(t,ii) *1000))*12/44;

file sol_state;
put sol_state;
sol_state.nw = 12;
sol_state.nr = 2;
sol_state.nz = 1e-15;
loop(t$(ord(t)<card(t)),
    put t.tl:4:0;
    loop(ii, put K.l(t,ii):11:3; );
    put MAT.l(t):11:3;
    put MU.l(t):11:3;
    put ML.l(t):11:3;
    loop(r, put TATM.l(t,r):11:3; );
    put TOCEAN.l(t):11:3;
    loop(ii, put AL.l(t,ii):11:3; );
    put /;
);


file sol_ctrl;
put sol_ctrl;
sol_ctrl.nw = 12;
sol_ctrl.nr = 2;
sol_ctrl.nz = 1e-15;
loop(t$(ord(t)<card(t)),
    put t.tl:4:0;
    loop(ii, put Cratio.l(t,ii):11:3; );
    loop(ii, put miu.l(t,ii):18:10; );
    put /;
);

file sol_other;
put sol_other;
sol_other.nw = 12;
sol_other.nr = 2;
sol_other.nz = 1e-15;
loop(t$(ord(t)<card(t)),
    put t.tl:4:0;

    loop(ii, put E.l(t,ii):11:3; );
    put F.l(t):11:3;
    loop(ii, put Y.l(t,ii):11:3; );
    loop(ii, put netY.l(t,ii):11:3; );
    loop(ii, put opt_SCC(t,ii):11:3; );
    loop(ii, put opt_carbonTax(t,ii):11:3; );
    put /;
);


**************************
parameter EnorthFixed(t)
    EtropicFixed(t)
;

positive VARIABLES
 Knorth(t)            Capital stock trillions US dollars
 MATnorth(t)          Carbon concentration in atmosphere GtC
 MUnorth(t)           Carbon concentration in shallow oceans Gtc
 MLnorth(t)
 Crationorth(t)            Consumption trillions US dollars
 miunorth(t)          Emission control rate GHGs
 Ynorth(t)
 ALnorth(t)
 netYnorth(t)
 mitiCostnorth(t)
 
 Ktropic(t)            Capital stock trillions US dollars
 MATtropic(t)          Carbon concentration in atmosphere GtC
 MUtropic(t)           Carbon concentration in shallow oceans Gtc
 MLtropic(t)
 Cratiotropic(t)            Consumption trillions US dollars
 miutropic(t)          Emission control rate GHGs
 Ytropic(t)
 ALtropic(t)
 netYtropic(t)
 mitiCosttropic(t)
;

variables
 Fnorth(t) 
 Enorth(t)            CO2-equivalent emissions GtC
 TATMnorth(t,r)         Temperature of atmosphere in degrees C
 TOCEANnorth(t)       Temperature of lower oceans degrees C
 UTILITYnorth
 
 Ftropic(t)
 Etropic(t)            CO2-equivalent emissions GtC
 TATMtropic(t,r)         Temperature of atmosphere in degrees C
 TOCEANtropic(t)       Temperature of lower oceans degrees C
 UTILITYtropic
;

*****************************************************************

EQUATIONS

 KKnorth(t)            Capital balance equation
 MMATnorth(t)          Atmospheric concentration equation
 MMUnorth(t)           Shallow ocean concentration
 MMLnorth(t)
 TATMEQ_Tropic_north(t)        Temperature-climate equation for atmosphere
 TATMEQ_North_north(t)        Temperature-climate equation for atmosphere
 TATMEQ_South_north(t)        Temperature-climate equation for atmosphere
 TOCEANEQnorth(t)      Temperature-climate equation for lower oceans
 TFPeqnorth(t)
 ForcingEqnorth(t)
 EEnorth(t)            Emissions equation
 MitigationCostFunnorth(t)
 Outputnorth(t)
 NetOutputnorth(t)
 UTILnorth             Objective function

 MMATtropic(t)          Atmospheric concentration equation
 KKtropic(t)            Capital balance equation 
 MMUtropic(t)           Shallow ocean concentration
 MMLtropic(t)
 TATMEQ_Tropic_tropic(t)        Temperature-climate equation for atmosphere
 TATMEQ_North_tropic(t)        Temperature-climate equation for atmosphere
 TATMEQ_South_tropic(t)        Temperature-climate equation for atmosphere
 TOCEANEQtropic(t)      Temperature-climate equation for lower oceans
 TFPeqtropic(t)
 ForcingEqtropic(t)
 EEtropic(t)            Emissions equation
 MitigationCostFuntropic(t)
 Outputtropic(t)
 NetOutputtropic(t)
 UTILtropic             Objective function
;

** Equations of the model

KKnorth(t)$(ord(t)<card(t))..       Knorth(t+1) =l= (1-DK)*Knorth(t) + (1-Crationorth(t))*netYnorth(t);
MMATnorth(t)$(ord(t)<card(t))..     MATnorth(t+1)  =e= MATnorth(t)*b11+MUnorth(t)*b21 + Enorth(t) + EtropicFixed(t) + ETREE(t);
MMUnorth(t)$(ord(t)<card(t))..      MUnorth(t+1)   =e= MATnorth(t)*b12+MUnorth(t)*b22+MLnorth(t)*b32;
MMLnorth(t)$(ord(t)<card(t))..      MLnorth(t+1)   =e= MUnorth(t)*b23+MLnorth(t)*b33;

TATMEQ_North_north(t)$(ord(t)<card(t))..     	TATMnorth(t+1,'North') =e= (1-B)*TATMnorth(t,'North') - 
						xi2*(TATMnorth(t,'North')-TOCEANnorth(t)) -
						xi4*(TATMnorth(t,'North')-TATMnorth(t,'Tropic')) + 
						(xi1+xi7)*Fnorth(t);
TATMEQ_Tropic_north(t)$(ord(t)<card(t))..     	TATMnorth(t+1,'Tropic') =e= (1-B)*TATMnorth(t,'Tropic') - 
						xi2*(TATMnorth(t,'Tropic')-TOCEANnorth(t)) - 
						xi4/2*(TATMnorth(t,'Tropic')-TATMnorth(t,'North')) - 
                        xi5/2*(TATMnorth(t,'Tropic')-TATMnorth(t,'South')) +
						(xi1+xi8)*Fnorth(t);
TATMEQ_South_north(t)$(ord(t)<card(t))..     	TATMnorth(t+1,'South') =e= (1-B)*TATMnorth(t,'South') - 
						xi2*(TATMnorth(t,'South')-TOCEANnorth(t)) -
						xi5*(TATMnorth(t,'South')-TATMnorth(t,'Tropic')) + 
						xi1*Fnorth(t);

TOCEANEQnorth(t)$(ord(t)<card(t))..   	TOCEANnorth(t+1) =e= TOCEANnorth(t) + sum(r, xi3*(TATMnorth(t,r)-TOCEANnorth(t))) +
                        xi3*(TATMnorth(t,'Tropic')-TOCEANnorth(t));


ForcingEqnorth(t)$(ord(t)<card(t))..   	Fnorth(t) =e= FCO22X*(log(MATnorth(t)/mateq)/log(2))+FORCOTH(t);

EEnorth(t)$(ord(t)<card(t))..          	Enorth(t) =e= sigma(t,'North') * (1-miunorth(t))*Ynorth(t);

TFPeqnorth(t)$(ord(t)<card(t))..		ALnorth(t) =e= AL_noDam(t,'North') / ( 1+sum(s$(ord(s)<=ord(t) and ord(s)>=ord(t)-nlag),
                                            (ga1('North')*(TATMnorth(t,'North')-TATM0('North'))+ga2('North')*sqr(TATMnorth(t,'North')-TATM0('North')))*delaT('North')**(ord(t)-ord(s))) );

Outputnorth(t)$(ord(t)<card(t))..		Ynorth(t) =e= (ALnorth(t)*Pop(t,'North')**(1-GAMA)*Knorth(t)**GAMA);
MitigationCostFunnorth(t)$(ord(t)<card(t)).. mitiCostnorth(t) =e= cost1(t,'North')*(miunorth(t)**EXPcost2);
NetOutputnorth(t)$(ord(t)<card(t))..	netYnorth(t) =e= (1-mitiCostnorth(t))*Ynorth(t);

UTILnorth..             UTILITYnorth =e= sum( (t)$(ord(t)<card(t)), (Crationorth(t)*netYnorth(t)/Pop(t,'North'))**elasmu1/elasmu1 * Pop(t,'North') * beta**(ord(t)-1) );


KKtropic(t)$(ord(t)<card(t))..       Ktropic(t+1) =l= (1-DK)*Ktropic(t) + (1-Cratiotropic(t))*netYtropic(t);
MMATtropic(t)$(ord(t)<card(t))..     MATtropic(t+1)  =e= MATtropic(t)*b11+MUtropic(t)*b21 + Etropic(t) + EnorthFixed(t) + ETREE(t);
MMUtropic(t)$(ord(t)<card(t))..      MUtropic(t+1)   =e= MATtropic(t)*b12+MUtropic(t)*b22+MLtropic(t)*b32;
MMLtropic(t)$(ord(t)<card(t))..      MLtropic(t+1)   =e= MUtropic(t)*b23+MLtropic(t)*b33;

TATMEQ_North_tropic(t)$(ord(t)<card(t))..     	TATMtropic(t+1,'North') =e= (1-B)*TATMtropic(t,'North') - 
						xi2*(TATMtropic(t,'North')-TOCEANtropic(t)) -
						xi4*(TATMtropic(t,'North')-TATMtropic(t,'Tropic')) + 
						(xi1+xi7)*Ftropic(t);
TATMEQ_Tropic_tropic(t)$(ord(t)<card(t))..     	TATMtropic(t+1,'Tropic') =e= (1-B)*TATMtropic(t,'Tropic') - 
						xi2*(TATMtropic(t,'Tropic')-TOCEANtropic(t)) - 
						xi4/2*(TATMtropic(t,'Tropic')-TATMtropic(t,'North')) - 
                        xi5/2*(TATMtropic(t,'Tropic')-TATMtropic(t,'South')) +
						(xi1+xi8)*Ftropic(t);
TATMEQ_South_tropic(t)$(ord(t)<card(t))..     	TATMtropic(t+1,'South') =e= (1-B)*TATMtropic(t,'South') - 
						xi2*(TATMtropic(t,'South')-TOCEANtropic(t)) -
						xi5*(TATMtropic(t,'South')-TATMtropic(t,'Tropic')) + 
						xi1*Ftropic(t);

TOCEANEQtropic(t)$(ord(t)<card(t))..   	TOCEANtropic(t+1) =e= TOCEANtropic(t) + sum(r, xi3*(TATMtropic(t,r)-TOCEANtropic(t))) +
                        xi3*(TATMtropic(t,'Tropic')-TOCEANtropic(t));


ForcingEqtropic(t)$(ord(t)<card(t))..   Ftropic(t) =e= FCO22X*(log(MATtropic(t)/mateq)/log(2))+FORCOTH(t);

EEtropic(t)$(ord(t)<card(t))..          Etropic(t) =e= sigma(t,'Tropic') * (1-miutropic(t))*Ytropic(t);

TFPeqtropic(t)$(ord(t)<card(t))..		ALtropic(t) =e= AL_noDam(t,'Tropic') / ( 1+sum(s$(ord(s)<=ord(t) and ord(s)>=ord(t)-nlag),
                                            (ga1('Tropic')*(TATMtropic(t,'Tropic')-TATM0('Tropic'))+ga2('Tropic')*sqr(TATMtropic(t,'Tropic')-TATM0('Tropic')))*delaT('Tropic')**(ord(t)-ord(s))) );

Outputtropic(t)$(ord(t)<card(t))..		Ytropic(t) =e= (ALtropic(t)*Pop(t,'Tropic')**(1-GAMA)*Ktropic(t)**GAMA);
MitigationCostFuntropic(t)$(ord(t)<card(t)).. mitiCosttropic(t) =e= cost1(t,'Tropic')*(miutropic(t)**EXPcost2);
NetOutputtropic(t)$(ord(t)<card(t))..	netYtropic(t) =e= (1-mitiCosttropic(t))*Ytropic(t);

UTILtropic..             UTILITYtropic =e= sum( (t)$(ord(t)<card(t)), (Cratiotropic(t)*netYtropic(t)/Pop(t,'Tropic'))**elasmu1/elasmu1 * Pop(t,'Tropic') * beta**(ord(t)-1) );

*******************************
**  Upper and Lower Bounds

Knorth.lo(t) = 0.01;
MATnorth.lo(t) = 10;
MUnorth.lo(t) = 1;
Ynorth.lo(t) = 0.01;
netYnorth.lo(t) = 0.01;

Crationorth.lo(t) = 0.1;
Crationorth.up(t) = 0.9;
miunorth.lo(t) = 0.001;
miunorth.up(t) = LIMMIU;
miunorth.up(t)$(ord(t)<150) = 1;
Enorth.lo(t)$(ord(t)<150) = 0;
TATMnorth.lo(t,r)$(ord(t)<150) = 0;

* Initial conditions
Knorth.FX(tfirst) = k0('North');
MATnorth.FX(tfirst) = mat0;
MUnorth.FX(tfirst) = mu0;
MLnorth.fx(tfirst) = ml0;
TATMnorth.FX(tfirst,r) = tatm0(r);
TOCEANnorth.FX(tfirst) = tocean0;
ALnorth.fx(tfirst) = A0('North');
    
Ktropic.lo(t) = 0.01;
MATtropic.lo(t) = 10;
MUtropic.lo(t) = 1;
Ytropic.lo(t) = 0.01;
netYtropic.lo(t) = 0.01;

Cratiotropic.lo(t) = 0.01;
Cratiotropic.up(t) = 1;
miutropic.lo(t) = 0.001;
miutropic.up(t) = LIMMIU;
miutropic.up(t)$(ord(t)<150) = 1;

* Initial conditions
Ktropic.FX(tfirst) = k0('Tropic');
MATtropic.FX(tfirst) = mat0;
MUtropic.FX(tfirst) = mu0;
MLtropic.fx(tfirst) = ml0;
TATMtropic.FX(tfirst,r) = tatm0(r);
TOCEANtropic.FX(tfirst) = tocean0;
ALtropic.fx(tfirst) = A0('Tropic');

model SocialProbNorth / KKnorth, MMATnorth, MMUnorth, MMLnorth, TATMEQ_South_north, TATMEQ_Tropic_north,
  TATMEQ_North_north, TOCEANEQnorth, TFPeqnorth, ForcingEqnorth, MitigationCostFunnorth, 
  EEnorth, Outputnorth, NetOutputnorth, Utilnorth /;
  
model SocialProbTropic / KKtropic, MMATtropic, MMUtropic, MMLtropic, TATMEQ_South_tropic, TATMEQ_Tropic_tropic,
  TATMEQ_North_tropic, TOCEANEQtropic, TFPeqtropic, ForcingEqtropic, MitigationCostFuntropic, 
  EEtropic, Outputtropic, NetOutputtropic, Utiltropic /;



Knorth.l(t) = K.l(t,'North');
MATnorth.l(t) = MAT.l(t);
MUnorth.l(t) = MU.l(t);
MLnorth.l(t) = ML.l(t);
TATMnorth.l(t,r) = TATM.l(t,r);
TOCEANnorth.l(t) = TOCEAN.l(t);
ALnorth.l(t) = AL.l(t,'North');

Ynorth.l(t) = Y.l(t,'North');
netYnorth.l(t) = netY.l(t,'North');
Enorth.l(t) = E.l(t,'North');
Fnorth.l(t) = F.l(t);
mitiCostnorth.l(t) = mitiCost.l(t,'North');

Crationorth.l(t) = Cratio.l(t,'North');
miunorth.l(t) = miu.l(t,'North');


Ktropic.l(t) = K.l(t,'Tropic');
MATtropic.l(t) = MAT.l(t);
MUtropic.l(t) = MU.l(t);
MLtropic.l(t) = ML.l(t);
TATMtropic.l(t,r) = TATM.l(t,r);
TOCEANtropic.l(t) = TOCEAN.l(t);
ALtropic.l(t) = AL.l(t,'Tropic');

Ytropic.l(t) = Y.l(t,'Tropic');
netYtropic.l(t) = netY.l(t,'Tropic');
Etropic.l(t) = E.l(t,'Tropic');
Ftropic.l(t) = F.l(t);
mitiCosttropic.l(t) = mitiCost.l(t,'Tropic');

Cratiotropic.l(t) = Cratio.l(t,'Tropic');
miutropic.l(t) = miu.l(t,'Tropic');


parameter iter 
    maxiternum / 1000 /
    omega / 0.5 /
    oldE(t,ii)
    currE(t,ii)
    err
;

currE(t,ii) = E.l(t,ii);

err = 1; 
iter = 0;
option solprint = off;
Repeat(        
    EnorthFixed(t) = currE(t,'North');    
     
    nnn = 0;
    nsucc = 0;
    Repeat(
        nnn = nnn+1;
        solve SocialProbTropic maximizing UTILITYtropic using nlp ;
        if((SocialProbTropic.MODELSTAT eq 2 and SocialProbTropic.SOLVESTAT eq 1),
          nsucc = nsucc+1;
        else
          nsucc = 0;
        );
    until (nnn>=nnnmax or nsucc>=nsuccmax) );
    
    ABORT$(nsucc < nsuccmax) "Model not normally completed";
    
    EtropicFixed(t) = currE(t,'Tropic');
    display EtropicFixed;
    
    nnn = 0;
    nsucc = 0;    
    Repeat(
        nnn = nnn+1;
        solve SocialProbNorth maximizing UTILITYnorth using nlp ;
        if((SocialProbNorth.MODELSTAT eq 2 and SocialProbNorth.SOLVESTAT eq 1),
          nsucc = nsucc+1;
        else
          nsucc = 0;
        );
    until (nnn>=nnnmax or nsucc>=nsuccmax) );
    
    ABORT$(nsucc < nsuccmax) "Model not normally completed";

    oldE(t,ii) = currE(t,ii);
    currE(t,'North') = (1-omega)*Enorth.l(t) + omega*oldE(t,'North');
    currE(t,'Tropic') = (1-omega)*Etropic.l(t) + omega*oldE(t,'Tropic');
    
    iter = iter+1;
    err = smax((t,ii)$(ord(t)<card(t)), abs(currE(t,ii)-oldE(t,ii))/(1+abs(oldE(t,ii))));
    display iter, err;
until(iter>=maxiternum or err<0.000001) );
    
ABORT$(err>0.000001) "No convergence!!!";

opt_SCC(t,'Tropic')$(ord(t)<card(t))=-1*MMATtropic.m(t)*1000/(KKtropic.m(t)+0.00000000001)*12/44;
opt_carbonTax(t,'Tropic')$(ord(t)<card(t))=(cost1(t,'Tropic')*(expcost2*miutropic.l(t)**(expcost2-1)/sigma(t,'Tropic') *1000))*12/44;

opt_SCC(t,'North')$(ord(t)<card(t))=-1*MMATnorth.m(t)*1000/(KKnorth.m(t)+0.00000000001)*12/44;
opt_carbonTax(t,'North')$(ord(t)<card(t))=(cost1(t,'North')*(expcost2*miunorth.l(t)**(expcost2-1)/sigma(t,'North') *1000))*12/44;


file sol_state_NE;
put sol_state_NE;
sol_state_NE.nw = 12;
sol_state_NE.nr = 2;
sol_state_NE.nz = 1e-15;
loop(t$(ord(t)<card(t)),
    put t.tl:4:0;
    put Knorth.l(t):11:3;
    put Ktropic.l(t):11:3;
    put MATnorth.l(t):11:3;
    put MUnorth.l(t):11:3;
    put MLnorth.l(t):11:3;
    loop(r, put TATMnorth.l(t,r):11:3; );
    put TOCEANnorth.l(t):11:3;
    put ALnorth.l(t):11:3;
    put ALtropic.l(t):11:3; 
    put /;
);


file sol_ctrl_NE;
put sol_ctrl_NE;
sol_ctrl_NE.nw = 12;
sol_ctrl_NE.nr = 2;
sol_ctrl_NE.nz = 1e-15;
loop(t$(ord(t)<card(t)),
    put t.tl:4:0;
    put Crationorth.l(t):11:3;
    put Cratiotropic.l(t):11:3;
    put miunorth.l(t):18:10;
    put miutropic.l(t):18:10; 
    put /;
);

file sol_other_NE;
put sol_other_NE;
sol_other_NE.nw = 12;
sol_other_NE.nr = 2;
sol_other_NE.nz = 1e-15;
loop(t$(ord(t)<card(t)),
    put t.tl:4:0;

    put Enorth.l(t):11:3;
    put Etropic.l(t):11:3;
    put Fnorth.l(t):11:3;
    put Ynorth.l(t):11:3;
    put Ytropic.l(t):11:3;
    put netYnorth.l(t):11:3; 
    put netYtropic.l(t):11:3; 
    loop(ii, put opt_SCC(t,ii):11:3; );
    loop(ii, put opt_carbonTax(t,ii):11:3; );
    put /;
);

